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MOLECULAR  DYNAMICAL  CALCULATIONS  OF  THE  THERMAL  DIFFUSIVITY  OF  A 
PERFECT  LATTICE* 

R.  A.  MacDonald  and  D.  H.  Tsai 
National  Bureau  of  Standards 
Washington,  D.C.  20234 


We  have  used  the  method  of  molecular  dynamics  to  make 
a detailed  study  of  thermal  dlffuslvlty  In  a perfect 
monatomic  lattice.  The  interatomic  potential  Is  that 
appropriate  for  iron.  We  limit  the  atomic  motions  to 
two  dimensions  In  order  to  shorten  the  computation. 

We  maintain  one  end  of  the  lattice  at  a given  kinetic 
temperature  and  obtain  the  temperature  profile  In  the 
lattice  as  a function  of  time.  The  total  energy  added 
to  the  system  is  recorded.  We  fit  diffusive  curves  to 
the  temperature  profiles  and  thus  obtain  the  thermal 
dlffuslvlty  of  the  lattice.  Its  value  is  4xl0“6  m2sec~1 
at  a mean  lattice  temperature  of  75K. 


1.  INTRODUCTION 

We  have  used  the  method  of  molecular  dynamics  to  study  the 
transport  of  thermal  energy  in  a perfect  monatomic  lattice.  The 
complexity  and  immensity  of  such  calculations  make  it  desirable  to 
relate  our  results  to  physical  reality  in  as  many  ways  as  possible. 
To  this  end,  we  obtained  a rough  estimate  of  the  lattice  thermal 


*This  work  was  performed  with  partial  support  from  the  U.S.  Energy 
Research  and  Development  Administration  and  from  the  U.S.  Army 
Research  Office,  Durham,  N.C. 
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conductivity  from  an  earlier  calculation  of  heat  pulae  propagation 
in  a three-dimensional  lattice^1).  This  led  to  the  present,  more 
detailed  study  of  thermal  diffuaivity  in  a lattice.  We  have  limit- 
ed the  atomic  motions  to  two  dimensions  in  order  to  shorten  the 
computation  time,  otherwise,  this  restriction  is  insignificant. 

We  shall  first  describe  our  model,  then  in  section  3 we  shall 
treat  the  problem  using  conventional  theory  of  heat  transfer  by 
diffusion  in  a continuum.  In  section  4 we  present  the  results  of 
our  molecular  dynamical  calculations  and  compare  them  with  the 
continuum  theory.  In  the  final  section  we  discuss  our  results. 


2.  MODEL 


Our  two-dimensional  system  is  made  up  of  bcc  unit  cells  forming 
a ribbon-like  filament  25  units  wide,  250  units  long,  and  1 unit 
thick.  This  ribbon  is  placed  in  Cartesian  coordinates  with  x de- 
noting the  width,  z the  length  and  y the  thickness.  Two-dimension- 
al motion  is  obtained  by  limiting  the  motion  in  the  y direction  to 
zero.  Mirror  boundary  conditions  are  Imposed  on  the  transverse 
boundaries  in  the  ± x directions  and  at  the  heated  end,  z - 0.  The 
interatomic  potential  Is  that  constructed  by  Chang^2)  for  a-iron. 
Initially,  this  lattice  is  in  thermal  equilibrium  at  a temperature 
T^.  Then,  the  first  ten  lattice  (x-y)  planes  are  heated  quickly  to 
an  average  temperature  T^  and  maintained  at  that  level  for  the  rest 
of  the  calculations.  The  amount  of  energy  added  to  these  ten  planes 
is  recorded.  The  classical  equations  of  motion  for  the  lattice 
atoms  are  solved  by  numerical  integration,  yielding  the  position 
and  velocity  of  each  atom.  The  average  kinetic  energy,  potential 
energy,  stress  components  and  density  of  lattice  planes  are  obtain- 
ed as  functions  of  time.  These  results  show  the  transport  of 
energy  into  the  lattice. 

We  use  the  rigid  mirror  boundary  condition  in  the  transverse 
direction  to  prevent  lateral  displacement  and  distortion  of  the 
filament.  This  was  necessary  because  we  found  that  when  periodic 
boundary  conditions  were  used,  as  in  earlier  calculations,  the 
filament  would  drift  laterally  from  the  equilibrium  position,  due 
to  accusmlation  of  small  errors  in  the  position  and  velocity  of  the 
atoms,  and  perhaps  also  dua  to  the  procedure  of  numerical  integra- 
tion. Since  this  motion  was  not  uniform  along  the  length  of  the 
filament,  and  since  the  undisturbed  end  was  constrained  to  stay  in 
place,  the  filament  became  laterally  distorted  and  strain  energy 
developed.  Over  the  duration  of  a run  (400  time  units)  the  strain 
energy  was  such  that  a two  fold  Increase  in  the  total  energy  was 
observed.  This  strain  energy  is  eliminated  when  the  mirror  bound- 
ary condition  is  applied. 
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3.  CONTINUUM  THEORY 

Since  the  lattice  is  at  a moderate  initial  temperature  (-0.1  0 , 
where  Q - Debye  temperature)  the  mechanism  of  heat  transfer  can  beD 
expected  to  be  diffusive.  It  is  therefore  appropriate  for  us  to 
test  our  results  against  those  of  the  continuum  theory  of  thermal 

diffusion^  - We  nrnreed  as  fnllnua* 


We  proceed  as  follows: 


When  a time-dependent  boundary  temperature  T,  (t)  is  imposed  on 
a system  initially  at  a uniform  temperature  T.,  "the  temperature 
profile  is  given  by  (see  ref.  3 pp.  273-4) 

T(z,t)  - TA  erf(u)  + (2 //»)  / Tfc(T-z2/4au2)exp(-u2)du  (1) 

•'ll 

where  u2  - z2/4a(r-t),  erf(u)  * (2 l^v)J  exp(-t2)dt, 

a is  the  thermal  diffusivity  and  z is  the  distance  in  the  direction 
of  heat  flow.  Although  we  maintain  the  average  value  of  the  kinetic 
energy  in  the  first  ten  planes  of  the  lattice  at  a constant  value, 
the  local  kinetic  tenqjerature  of  the  tenth  (boundary)  plane  varies 
considerably  on  an  atomic  scale.  For  the  diffusive  calculation  we 
need  the  average  temperature  of  the  boundary  plane.  In  the  present 
calculation,  this  may  be  approximated  by  three  linear  segments  (see 
7ig.  1): 

T.(t)  " T0+  ciT»  0 < 


VT>  * Tl+  C2t  1 

Tb(T)  ■ T2» 


Tj<  T $ T2 


t2<  T, 


Substitution  into  Eq.  1 yields  the  following  expression  for  the 
temperature  profile: 

T(z,t)  - Tj+  [Tq-  Tt+  CjT]F(uo)  - [Tq-  T^  (Cj-  C2)t]F(Ui) 

f (T2-  Tr  C2t]F(u2) 

+ 2 C.  t [u  F(u2)  - u exp(-u2)//*] 
loo  o o 

- 2 (C j—  C2)  (t-Tj)  (u2F(uj)  - Ujexp(-u2)//ir] 

- 2C2(t-t2)[u|f(u2)  - u2exp  (-u2)//i»] . C 

where  F(u)-  1 - erf(u),  uQ  - z/2*''ar', 

Uj  - z/2/o(t-Tj) , u2  ■ x/2/o(t-t2). 
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T(z,t)  can  be  fitted  to  our  computed  temperature  profiles  at  one 
value  of  t.  Since  the  diffusive  profile  is  given  in  terms  of 
u(“z/2*^or) , the  distance  scale  is  arbitrary  when  a is  unknown. 

This  arbitrariness  is  removed  by  requiring  that  the  area  under  the 
diffusive  curve  be  equal  to  the  total  energy  a^ded  to  the  system  up 
to  the  time  t.  With  the  value  of  a thus  obtained,  the  diffusive 
curves  at  other  times  are  obtained  by  appropriately  scaling  the 
distance.  The  area  under  these  other  diffusive  curves  must  also  be 
in  agreement  with  the  total  energy  added  to  the  system  at  their 
respective  times.  If  this  is  so,  it  can  be  concluded  that  the  heat 
transfer  in  our  system  is  indeed  diffusive. 


4.  RESULTS 


In  Fig.  1,  we  show  the  boundary  temperature  as  a function  of 
time,  T,  (t).  Even  though  we  have  averaged  the  kinetic  temperature 
of  planes  6 to  15  to  obtain  T^(t),  there  is  still  considerable 
fluctuation  in  the  boundary  temperature  and  we  approximate  it  by 
the  dashed  curve  C,  as  given  in  Eq.  2.  The  constants  in  these 
equations  have  the  following  values: 

Tq  - 59  K,  Tx  - 102  K,  Tj  - 115  K, 

t - 20,  t2  - 208,  Cj  - 2.214  K and  C2  - 0.0628  K. 

The  time  t is  measured  in  units  of  lattice  plane  spacing  d divided 
by  longitudinal  sound  velocity  (equal  to  0.264xl0“^8)  • The  initial 
temperature  of  the  lattice,  T.,  is  44  K.  These  values  are  used  in 
Eq.  3 to  calculate  the  diffusive  temperature  profiles. 

In  Fig.  2 we  show  the  kinetic  temperature,  density  and  longitu- 
dinal stress  profiles  in  the  lattice  at  time  t - 300.  Each  point 
is  the  average  of  10  neighboring  planes.  The  kinetic  temperature 
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Fig.l.  Boundary  temperature  (Kelvin)  as  a function  of  time  t (unit 
of  t - 0.264xl0“13s).  Tfe  is  the  average  kinetic  temperature 
of  planes  6 to  15.  TheDboundary  is  at  plane  10.  C labels 
the  approximate  value  of  boundary  temperature  used  in  the 
continuum  theory.  G denotea  a gap  in  computer  output. 
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Fig. 2.  Kinetic  temperature  (Kelvin),  relative  density  and  longitudi- 
nal stress  (Pascal)  versus  lattice  plane  number.  Each  point 
is  the  average  over  10  neighboring  planes.  For  the  density 
p/p  and  stress,  S , t - 300;  for  the  kinetic  temperature, 
t -°298  (A)  and  t 5*303  (9).  Continuum  theory  diffusive  pro- 
file, C,  is  at  t - 300  ( ).  T£  « 44  K is  the  uniform 

initial  temperature  of  the  lattice. 

of  a plane  is  defined  as  m<V2>/2k^  where  m is  the  mass  of  an  atom, 

V its  velocity,  <•••>  indicates  the  average  over  all  atoms  in  the 
plane,  and  is  Boltsmann's  constant.  We  plot  two  temperature 
profiles  at  slightly  different  times,  t - 298  and  t - 303,  to  show 
how  much  local  fluctuation  there  is  in  the  temperature.  The  diffu- 
sive curve,  C,  is  calculated  frqp  $q.  3 at  t ■ 300.  The  prominent 
pulses  in  the  stress  and  density  profiles  near  lattice  plane  300 
are  due  to  the  stress  and  strain  in  the  heated  planes  (and  their 
image  across  the  mirror)  generated  during  the  initial  rapid  heating. 
These  pulses  propagate  with  the  longitudinal  velocity  of  sound,  as 
may  be  determined  from  their  trajectories  in  a series  of  profiles. 
Although  there  is  a slight  drop  in  pressure  near  the  heated 
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Fig. 3.  Energy  in  lattice  (eV)  as  a function  of  time  t.  D is  the 

energy  added  to  the  heated  planes,  per  filament  of  unit  cell 
cross  section.  G denotes  a gap  in  computer  output.  C is 
energy  added  to  the  filament  as  calculated  by  continuum 
theory  of  diffusion  with  a - 4xl0~6  m2  s"1.  E is  total 
energy  in  filament  minus  its  initial  energy  content. 


boundary  at  this  particular  time  (t  - 300),  our  data  show  that  both 
the  time  average  and  the  spatial  average  of  the  stress  have  a 
constant  value  of  1. 84x10*3  Pascal  (1.84  kllobar)  behind  the  initial 
pulse.  This  indicates  that,  as  expected,  the  stresses  equilibrate 
with  the  velocity  of  sound.  The  density  profile  shows  that  the 
density  near  the  heated  end  is  definitely  lower  than  the  average 
density.  This  is  clearly  due  to  the  higher  temperature  at  the 
heated  end  and  is,  in  our  view,  responsible  for  the  observed  fall 
of  the  computed  temperature  profile  below  the  diffusive  curve  C 
near  the  heated  end  of  the  lattice.  These  results  therefore  show 
quite  clearly  the  thermoelastic  coupling  in  the  lattice,  a feature 
that  is  usually  Ignored  in  the  continuum  theory  of  thermal  diffusion. 

In  Fig.  3,  curve  D is  the  energy  added  to  the  lattice  at  the 
heated  end  per  filament  of  wit  cell  cross-section  as  a function  of 
time.  We  calculate  the  diffusivlty  of  the  lattice  from  the  diffu- 
sive curve  C at  r • 300  in  Pig.  2,  by  requiring  that  the  area  under 
that  curve  be  equal  to  the  value  of  curve  D at  t • 300  in  Fig.  3, 
as  described  in  section  3.  We  obtain  the  following  diffusivlty, 

a " 4.0x10”*  n2  s-1, 
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or  conductivity 


k - 9.4  W m-1  K_1 


where  we  have  used  the  relation 


k **  apC 


akB/2d3. 


We  use  this  value  of  a to  obtain  diffusive  curves  at  other  times 
from  Eq.  3.  The  areas  under  these  curves  are  plotted  as  curve  C 
in  Fig.  3,  in  good  agreement  with  the  computed  curve  D. 

The  curve  E shows  the  total  energy  that  is  actually  in  a 
filament  of  the  lattice  minus  the  energy  that  was  in  the  lattice  at 
uniform  temperature  T,  before  heating.  The  difference  between 
curves  D and  E shows  the  lack  of  energy  conservation  in  the  system 
due  to  truncation  error  in  the  numerical  procedure.  By  the  end  of 
our  run,  this  error  amounts  to  about  3%  of  the  total  energy  of  the 
system. 


5.  DISCUSSION 

The  present  calculation  is  an  outgrowth  of  our  earlier  study 
of  the  propagation  of  a heat  pulse  in  a lattice.  One  of  our  major 
objectives  here  is  to  seek  additional  confirmation,  beyond  what  has 
been  established  in  previous  calculations,  that  our  model  is  indeed 
realistic.  Toward  this  end,  we  are  able  to  show  that  the  heat  flow 
in  our  model  is  diffusive,  as  we  expected.  There  is  no  experimental 
value  of  the  lattice  thermal  conductivity  of  pure  crystalline  iron 
with  which  to  compare  our  value  of  <,  but  for  the  alloy  Fe  (99.5%) 

Ni  (0.5%)  there  is  a recent  estimate  viz.,  < «*  28.2  W m_1  K_1  at 
75  r(4).  Considering  the  fact  that  the  effective  interatomic 
potential  for  iron  is  not  known  at  all  accurately,  we  feel  that  the 
agreement  obtained  here  is  really  quite  satisfactory.  In  addition, 
our  results  provide  a wealth  of  details  for  example,  we  see  that 
there  is  thermoelastic  coupling  in  the  process  of  heating,  and  there 
is  energy  non-conservation  due  to  the  unexpected  strain  energy  in 
the  lattice  caused  by  the  cumulative  errors  in  the  numerical  proce- 
dure. In  this  connection,  we  are  able  to  show  that  when  mirror 
boundary  conditions  are  Imposed  the  residual  numerical  error  is  in- 
deed small,  a few  percent.  There  is  also  evidence  of  second  sound, 
albeit  heavily  damped,  propagating  at  this  moderately  high  tempera- 
ture. However,  we  shall  not  pursue  this  intriguing  problem  here. 
Further  tests  of  the  model  may  be  devised.  For  example,  it  would  be 
useful  to  Investigate  the  effect  of  temperature  on  the  lattice 
thermal  conductivity.  If  the  result  should  also  prove  to  be  reason- 
able, it  would  then  be  of  interest  to  investigate  thermal  conductiv- 
ity under  conditions  not  easily  accessible  by  present-day  theory  or 
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experiment,  e.g.,  under  conditions  of  extreme  pressure.  Our  expe- 
rience indicates,  however,  that  the  calculation  of  thermal  diffu- 
sivity  by  the  method  of  molecular  dynamics  is  so  lengthy  as  to  be 
too  costly  and  too  inefficient.  If  the  transient  nature  of  the 
heat  flow  is  not  under  study,  then  we  believe  that  a model  for 
steady  heat  flow,  from  which  the  thermal  conductivity  may  be  obtain- 
ed directly (5),  would  be  considerably  more  economical  and  more 
convenient  to  investigate  than  the  diffusive  model  employed  here. 
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